#####################################################################################
#                                                                                   #
#    Scripts for breeding songbird point count data in New Jersey Pine Barrens      #
#                                                                                   #
#                           Created by: Philip M. Coppola                           #
#                                                                                   #
#####################################################################################

##### FORMATING DATA FROM ORIGINAL ENTRY SPREADSHEETS TO SINGLE CSV FOR ANALYSIS

# Install Packages
#install.packages(c("chron", "suncalc", "plyr", "stringr"))

# Load Packages
packagelist <- c("chron", "suncalc", "plyr", "stringr")
lapply(packagelist, require, character.only = TRUE)

# Set working directory - update as needed
setwd("E:/My Drive/Phil Coppola Dissertation Research files/All Data/Coppola NJPB Songbird/Outputs/")

# Import songbird count and site characterisitcs data and attach column names
songall <- as.data.frame(read.csv("../Combined/SONG_NJPB_ALL_MASTER_Distance-TR-Counts.csv", header = TRUE))
loc <- as.data.frame(read.csv("../Combined/SONG_NJPB_ALL_MASTER_PointLocations.csv", header = TRUE))
#veg <- as.data.frame(read.csv("../Combined/SONG_NJPB_ALL_MASTER_Vegetation.csv", header = TRUE))
attach(songall)

# Create 8 new columns combing M, F, and U for each combination of distance bin (2) and time interval (4); thus ignoring sex of indv.
songall$aa <- songall$M11+songall$F11+songall$U11  ;  songall$ab <- songall$M12+songall$F12+songall$U12
songall$ba <- songall$M21+songall$F21+songall$U21  ;  songall$bb <- songall$M22+songall$F22+songall$U22
songall$ca <- songall$M31+songall$F31+songall$U31  ;  songall$cb <- songall$M32+songall$F32+songall$U32
songall$da <- songall$M41+songall$F41+songall$U41  ;  songall$db <- songall$M42+songall$F42+songall$U42

# Create new dataframe, that removes sex columns; then confirm that these values are either 0 or 1
song <- songall[c(1:24,49:56)]
min(song[25:32])  ;  max(song[25:32])

# Add column with detection distance and interval combination (1-8)
# 1 is detection in first distance bin within first time interval, 2 is second distance bin in first time interval, etc.
get.first <- function(x) min(which(x!=0))
song$det <- apply(song[25:32], 1, get.first)

# Add column with time removal interval (1, 2, 3, or 4) and distance (1 or 2) derived from det above
song$int[song$det == 1] <- 1  ;  song$dist[song$det == 1] <- 1
song$int[song$det == 2] <- 1  ;  song$dist[song$det == 2] <- 2
song$int[song$det == 3] <- 2  ;  song$dist[song$det == 3] <- 1
song$int[song$det == 4] <- 2  ;  song$dist[song$det == 4] <- 2
song$int[song$det == 5] <- 3  ;  song$dist[song$det == 5] <- 1
song$int[song$det == 6] <- 3  ;  song$dist[song$det == 6] <- 2
song$int[song$det == 7] <- 4  ;  song$dist[song$det == 7] <- 1
song$int[song$det == 8] <- 4  ;  song$dist[song$det == 8] <- 2

# Create new dataframe with time removal interval (1, 2, 3, or 4) and distance (1 or 2); then confirm max is 4 and 2, respectively
song <- song[c(1:24,34:35)]
max(song[25])  ;  max(song[26])

# Format date in dataframe
song$date <- as.Date(song$date, "%m/%d/%Y")
date <- sort(unique(song$date, fromLast = TRUE))
length(date)  ;  sort(date)

# Merge lat long data with point count data, derive sunrise time of surevy morning based on location and date
song <- merge(song, loc, by.x = "rtpt", by.y = "rtpt", all.x = TRUE)
sun <- data.frame(date = as.Date(song$date), lat = as.numeric(song$lat), lon = as.numeric(song$long))
rise <- getSunlightTimes(data = sun, keep = c("sunrise"), tz = "America/New_York")
song$timestamp <- as.POSIXct(paste(song$date, song$stime, sep = " "), format="%Y-%m-%d %H:%M")
song$rise <- as.POSIXct(rise$sunrise, format = "%H:%M:%S")

# Caulcate column containing minutes after sunrise for each point count; trim dataframe to remove intermediate columns
song$minasun <- as.numeric(round(difftime(song$timestamp, song$rise, units = c("mins"))))
song$jdn <- format(song$date, "%j")
song$jdn <- as.numeric(song$jdn)
song_base <- song[c(1:3,27:29,14:18,34:35,19:26,4:13)]

# Combine site characteristic data with point count data
#song_veg <- merge(song, veg, by = c("rtpt", "year"), all.x = TRUE)

# Otionally write .csv containing formatted point count and habitat characteristic data
#write.csv(song_base, "song_base.csv", row.names=FALSE)


#######################

##### Seperating by species

# Print sorted list of all species deteced (4-letter AOU Code)
spp <- sort(unique(as.vector(song_base$sp))); length(spp); sort(spp)

# Create a list of data frames for each species
spp_list = split(song_base, f = song_base$sp)

# Optional: Write csv for each species detection
#mapply(write.csv, spp_list, file=paste0("../Outputs/Species/", names(spp_list), ".csv"))

# Summary and descriptive statistics
summary(song_base)
stime <- sort(unique(c(as.matrix(song_base$stime)), fromLast = TRUE)); length(stime); sort(stime)
etime <- sort(unique(c(as.matrix(song_base$etime)), fromLast = TRUE)); length(etime); sort(etime)
sky <- sort(unique(c(as.matrix(song_base$sky)), fromLast = TRUE)); length(sky); sort(sky)
wind <- sort(unique(c(as.matrix(song_base$wind)), fromLast = TRUE)); length(wind); sort(wind)
disturb <- sort(unique(c(as.matrix(song_base$disturb)), fromLast = TRUE)); length(disturb); sort(disturb)
temp <- sort(unique(c(as.matrix(song_base$temp)), fromLast = TRUE)); length(temp); sort(temp)
hum <- sort(unique(c(as.matrix(song_base$hum)), fromLast = TRUE)); length(hum); sort(hum)
sites <- sort(unique(c(as.matrix(song_base$rtpt)), fromLast = TRUE)); length(sites); sort(sites)
year <- sort(unique(c(as.matrix(song_base$year)), fromLast = TRUE)); length(year); sort(year)
obs <- sort(unique(c(as.matrix(song_base$obs)), fromLast = TRUE)); length(obs); sort(obs)
visit <- sort(unique(c(as.matrix(song_base$visit)), fromLast = TRUE)); length(visit); sort(visit)

